Growth and Structure of Stochastic Sequences 
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We introduce a class of stochastic integer sequences. In these sequences, every element is a 
sum of two previous elements, at least one of which is chosen randomly. The interplay between 
randomness and memory underlying these sequences leads to a wide variety of behaviors ranging 
from stretched exponential to log-normal to algebraic growth. Interestingly, the set of all possible 
sequence values has an intricate structure. 
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Integer sequences underly many problems in combina- 
torics, computer science, and physics, with new beautiful 
sequences continuing to emerge Q. Sequences are typi- 
cally deterministic. Meanwhile, stochastic sequences are 
just as ubiquitous, occurring in random processes such 
as the random walk. Stochastic sequences usually arise 
in very different contexts, and hence are rarely compared 
with their deterministic counterparts. In this article, we 
demonstrate how rich such a comparison can be. 

Consider the Fibonacci numbers, F n = F n ^\ + F n —2, 
that describe for example the number of leaves in plants 
and the number of ancestors of a drone J^, 0, . As every 
element depends on the previous two, a natural stochastic 
generalization is x n = x„_i ± x„_2, where addition and 
subtraction are chosen with equal probabilities || 0, [t]. ||] 
(similar sequences also describe one-dimensional disor- 
dered systems [p[ fL0|). The resulting sequences are in- 
triguing. While the sequences still grow exponentially, 
the ratio x„/x„_i approaches a stationary distribution 
that possesses singularities at all rational values 0, ^ . 

Inspired by this richness, we consider an alternative 
form of stochasticity, namely, one that does not require 
subtraction and therefore more similar in spirit to the 
original deterministic sequence. Relaxing the rule that 
every element depends only on the preceding two ele- 
ments, we arrive at the following additive stochastic rules 



(model I); 
(model II). 



(1) 



In model I, we take the preceding element x n -i and an- 
other one whose index q is randomly chosen between 
and n — 1. In model II, both indices p and q are chosen 
randomly, < p, q < n — 1. Without loss of generality, 
the first element is set to unity, xq = 1. Consequently 
x\ = 2, while the next elements are stochastic. The 
number of possible sequences increases as n! and n! 2 for 
models I and II, respectively. Rule I leads to monoton- 
ically increasing sequences; sequences generated by rule 
II increase only on average. 



The most basic characteristic, the average of the nth 
element, A be determined analytically. For 

model I, it satisfies the linear recursion relation 



n— 1 

A n = A n -i + -Y\A j . 



(2) 



Comparing this with the recursion for A n+ \, we 
eliminate the summation and obtain a Fibonacci- 
like recursion relation with n-dependent coefficients 
A n+1 — 2A n + A n _i = (n + l) _1 ^4„_i. We are primarily 
interested in the large-n behavior, and hence, we treat A 
and n as continuous variables. The above difference equa- 
tion reduces to the differential equation A" = (A — A')/n 
(here ' = d/dn). Using the WKB method [|ll|, we obtain 
the n-dependence 

- 1 / 4 exp(2V^)- 



an 



(3) 

The amplitude a ~ 0.1711 is determined numerically. We 
see that the long range memory leads to considerably 
slower stretched exponential growth compared with the 
exponentially growing Fibonacci numbers. 

Does the average characterize the growth of an actual 
sequence? If yes, this would imply that the normalized 
moments (x*) / '(x n ) k approach finite values asymptoti- 
cally. Figure 1 shows otherwise: the higher order mo- 
ments grow according to 



oc 



exp (PkVn) , 



(4) 



with (3k > kfti (for the lowest moments, we find 0k = 2, 
4.3, and 6.5). This so-called "multiscaling" indicates that 
a typical sequence may greatly depart from the average. 
Therefore, the average (|||) is insufficient to describe a 
typical sequence. These results are similar in spirit to the 
behavior of the random Fibonacci sequences where the 
typical growth is x oc exp(7n), with nontrivial Lyapunov 
exponent 7, while (x^) oc exp^n) with jk 7^ &7i- 

Interestingly, an individual realization grows slower 
than the average (see the inset to Fig. 1) 



x n oc exp (fiy/n ) , 



(5) 
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FIG. 1: The moments {x^) 1 ^ versus n for model I. The mo- 
ments were obtained from an average over 10 8 realizations. 
The inset compares the growth of an individual realization 
x n with the average A n . 



FIG. 2: The scaling distribution underlying P n (In x n ). Shown 
is the scaling function <fr(z) versus the scaling variable z = 
(\nx n — b n )/A n . The distributions were calculated from 10 7 
realizations. A Gaussian is also shown for reference. 



with the Lyapunov exponent f3 « 1.889. This coefficient 
was determined by studying the variable lnx„. As shown 
in Fig. 2, this variable is Gaussian distributed 



P„(lnx„) oc exp 



(lnxn - KY 
2A2 



(6) 



The average and the variance of lna;„ grow with n ac- 
cording to b n ~ fin 1 / 2 and A 2 ~ a 2 v}^ 2 oc b n , respec- 
tively. Eventually, the random variable y = hiXn/n 1 / 2 
becomes deterministic, y — > (3 as n — > oo. Similar behav- 
ior, including the Gaussian fluctuations and the relation 
between the variance and the average, is also found in 
one-dimensional localization problems (12, 13 



One can calculate the probability distribution P n (x) 
for extremal values of x n . The minimal value n + 1 is 
obtained by choosing q = at every step k — 1, . . . , n. 
Similarly, the maximal value 2™ is obtained by choosing 
q = k — 1 at every step. Hence, P n (n + 1) = P n (2 n ) — 
l/nl. Further extremal cases can be evaluated manually, 
and for example, P n (n + 2) = P n (3 • 2 n ~ 2 ) = 2=1. How- 
ever, these extremal probabilities do not elucidate the 
typical behavior (H). 

It is interesting to study the set of all possible sequence 
values. Let f2„ be the support of the probability distri- 
bution P n {x). For small n, we have SIq = {1}, fii = {2}, 
and furthermore, 

n 2 = {3,4} 

n 3 = {4, 5, 6,*, 8} 

fi 4 = {5,..., 10,*, 12, *,*,*, 16} 

ft 5 = {6,..., 18,*, 20, *,*,*, 24, *,...,*, 32}. (7) 

Determination of the sets f2„ requires enumeration of all 
n\ histories, and we computed them up to n — 15. The 



simplest feature is the set size T„, listed in table I. The set 
VL n always begins with a subsequence of B n consecutive 
integers. For n > 3, the sets fi„ have gaps, i.e., strings of 
missing elements denoted by * in Eq. (7j). The number 
of gaps G n is listed in table I. The three sequences r„, 
B n , and G n , all grow exponentially with n. For example, 
r„ oc A™ with A 1.78. We conclude that the sets 
£l n contain a number of nontrivial deterministic integer 
sequences including r„, B n , and G n . 

Remarkably, the sets f2„ have an intricate structure. 
For instance, for n = 8 the sequence of the gap lengths 
contains Gs = 18 elements as follows 



{1,1,1,1,1,3,1,1,3,3,7, 7,1,3,7,15,31,63}. 



(8) 



Generally, going in reverse direction (from 2™ to n + B n ) 
one observes a family of consecutive gaps of lengths 
2n-2 _ ^ 2"~ 3 — 1, . . . , 1 separated by single elements. 
Then, there is a 3-element sequence, followed by a second 
family of twin gaps, 2 n ~ 5 — 1, 2 n ~ 5 — 1, . . . , 1, 1. All these 
gaps are separated by single elements. Next, there is a 
5-sequence, followed by a family of triplet gaps, again sep- 
arated by single elements [this family has not yet formed 
for n = 8, Eq. (g)]. There is also a fourth family of 
quadruplet gaps with an intertwined pattern. The com- 
plexity of this gap-sequence structure increases rapidly, 
and eventually, gaps of even length appear. 

Naively, one may probe the probability distribution 
via a mean-field description that ignores the sequence 
history altogether. In this approximation, one obtains a 
recursive equation for the probability distribution 



Pn(x) 



(9) 



i=o y=i 



3 




Consequently, there are closed recursion relations for the 
moments. While consistent with the exact equation (|^), 
the emerging recursion relations for the higher moments 
are only approximate. Analysis of these equations results 
in ordinary scaling behavior, (x™) oc (x n ) k , contrary to 
Eq. (^). Therefore, strong correlations develop, correla- 
tions that affect the statistical characteristics. 

We now turn to model II. Here, the average A n satisfies 
a recursion relation similar to Eq. (H) , 



A,, 



n-1 
3=0 



(10) 



Simplifying this equation to A n = (l + n x ) A n _i and 
using the initial condition Aq = 1, we obtain 



A n = n + 1. 



(11) 



Numerical simulations confirm this linear growth. Fig- 
ure 3 shows that the normalized moments remain finite 
asymptotically as 



(12) 



In contrast with model I, the average properly charac- 
terizes higher order moments. Therefore, the probability 
distribution P n (x) admits the scaling form 



P n (x) ~ n <$>(z), Z = xn 



(13) 



in the asymptotic limit x, n — > oo with the scaling vari- 
able z = xrT 1 fixed (see Fig. 4). 

This scaling behavior enables quantitative characteri- 
zation of extremal statistics. Numerically, we find that 
the scaling function exhibits the following extremal be- 
haviors 



$(z) (X 



z z — > 0, 

exp (— z K ) z — > oo, 



(14) 



with k » 0.4 — 0.5 (see the inset to Fig. 4). The small- 
z behavior can be understood by considering the mini- 
mal value x n = 2. This occurs when the first element 



FIG. 3: The moments {x k n ) 1/k versus n for model II. The data 
represents an average over 10 8 realizations. 



xq = 1 is chosen twice and therefore, P„(2) = n~ 2 . 
Combining this with Eq. (p"3|), gives <I>(2rt _1 ) = n^ 1 , in 
agreement with the asymptotics Q(z) oc z in the z — > 
limit. The large- z behavior is more subtle as it depends 
on the entire sequence evolution. Contrary to the small 
argument behavior, analysis of the maximal sequences 
x n = 2™ does not elucidate the large argument tail. In- 
deed, such sequences occur with probability 1/n!, much 
smaller than the exponentially small probabilities domi- 
nating the large-z behavior. 

The ordinary scaling behavior indicates that mean- 
field theory may provide better insight in the case of 
model II. Ignoring the history by which sequences evolve 
yields the following recursion for the distribution 



n— 1 7i— 1 x— 1 
1=0 m=0 y=l 



y)- 



(15) 



Using Eq. (|13j ), and replacing summations by integration, 
equation (p.5[) reduces into the integral equation 



1 d£ f 1 drj 



Uo 1 



dz'<S> 



i] 



(16) 



The convolution structure suggests using the Laplace 
transform, and indeed, F(s) — J dze~ sz §{z) satisfies a 
simple equation 



F{s) 



1 -|2 



(17) 



The auxiliary function G(s) = ds' F(s') obeys the or- 
dinary differential equation dG/ds = (G/s) 2 from which 
G(s) = s/(l+cs) and then F(s) = (l + cs)~ 2 . The small- 
s behavior F(s) = 1 — s implies c — 1/2. Inverting the 
Laplace transform F(s) = (1 + s/2) -2 yields the scaling 
function 



$(z) = 4z exp (— 2z) 



(18) 



4 




z 



FIG. 4: The scaling distribution <f>(z) vs. z. The distributions 
were obtained from an average over 10 7 realizations. Shown 
also is the mean field theory (MFT) result @. The inset 
shows the large argument tail. 

In the small- z limit, mean-field theory is correct because 
memory is irrelevant for x <C n, (see Fig. 3). In contrast, 
for x 3> n, memory is important and the exponent k = 1 
is larger than the numerical value k m 0.4—0.5. Addition- 
ally, one may compare the prefactors characterizing the 
moments defined in Eq. (|l2|). Using /i^ = j dzz k &(z) 
gives ^ifc = (k + l)!2~ fc and in particular fi^ — 1, 3/2, 3, 
for k = 1, 2, 3. The corresponding numerical values are 
1, 1.84, 5.76, respectively. Although the history indepen- 
dent approximation is quantitatively inaccurate, it still 
provides useful insights for model II. 



We have seen that the sequence growth sensitively de- 
pends on the details of the model. In fact, the stretched 
exponential or algebraic growth can be tuned by varying 
the recurrence rules. For example, introducing a multi- 
plicative factor to model II, x n = c(x p + x q ), leads to 
the algebraic growth A n ~ n 2c_1 . Functionally differ- 
ent growth laws naturally emerge as well. If in model I, 
x n — x n -\ + x q , the memory range is < q < [bn] with 
< b < 1, one finds log- normal growth 

A n oc exp (C In 2 n) , (19) 

with C = [21n(l/6)p 1 . This growth law is slower than 
stretched exponential but faster than power law. 

In summary, we have introduced a class of stochas- 
tic integer sequences where each sequence element is a 
sum of two previous elements, at least one of which is 
randomly chosen. While the sequence may attain a vast 
range of possible values, the dynamics chooses a much 
narrower range of values. Depending on the governing 
rules, there is a wide spectrum of growth from algebraic 
to log-normal to stretched exponential. In model I, there 
are infinitely many relevant scales underlying the mo- 
ments. In contrast, for model II, there is a single scale 
and consequently, a mean-field approximation is qualita- 
tively correct. 

Generally, the phase space has an intricate structure. 
It contains alternating sequences of consecutive integers 
marking allowed and forbidden sequence values. The gap 
structure consists of increasingly complex patterns. In- 
teresting deterministic integer sequences such as the size 
of the phase space, the number of gaps, and the size of 
the first accessible sequence underlie this phase space. 
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